// whgiirs2grid.cpp : ¶¨Òå¿ØÖÆÌ¨Ó¦ÓÃ³ÌÐòµÄÈë¿Úµã¡£
//
 
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>

#include "gdal_priv.h"

#include "wGdalRaster.h"

using namespace std;

#define MOD1 1
#define MOD2 2
#define MOD_4B 3
#define MOD_4B_AQ 4
 

enum USETRACKTYPE{
	TRACK_ASC,
	TRACK_DES,
	TRACK_BOTH
};
 


enum USESIDETYPE{
	USE_SIDE_FROM_0_TO_ROWX,    //Ê¹ÓÃÉÏÃæ²¿·Ö
	USE_SIDE_FROM_ROWX_TO_YSIZE,//Ê¹ÓÃÏÂÃæ²¿·Ö
	USE_SIDE_FROM_UNKNOWN
};
 

float wmax( float a , float b )
{
	if( a > b ) return a ;
	else return b ;
}

float wmin( float a , float b )
{
	if( a < b ) return a ;
	else return b ;
}

int wmaxi( int a , int b )
{
	if( a > b ) return a ;
	else return b; 
}
int wmini( int a , int b )
{
	if( a < b ) return a;
	else return b ;
}


//¾­¶È×ªÍø¸ñ
int lon2grid( float lon , float reso , float lonLeft )
{
	return (lon - lonLeft ) / reso   ;
}

float grid2lon( int ix , float reso , float lonLeft )
{
	return ix * reso + lonLeft + reso/2.f ;
}

float lon2gridlon( float lon , float reso  )
{
	int ix = (lon + 180.f ) / reso   ;
	return ix * reso - 180.f + reso/2.f ;
}

//Î³¶È×ªÍø¸ñ
int lat2grid( float lat , float reso , float latTop )
{
	return (latTop - lat ) /reso   ;
}
float grid2lat( int iy , float reso , float latTop )
{
	return latTop - iy * reso - reso/2.f ;
}
float lat2gridlat( float lat , float reso )
{
	int i = (90.f - lat ) / reso  ;
	return 90.f - i * reso - reso/2.f ;
}




     
 

struct llPoint{
	double lon,lat;
	int x,y;
	vector<float> values;
	static int gridxsize ;
};
int llPoint::gridxsize = 0 ;

//¾­Î³¶Èµã
struct sllPoint
{
	float lon,lat ;
	inline sllPoint(float tlon,float tlat):lon(tlon),lat(tlat){};
	inline sllPoint():lon(0),lat(0){} ;
};


struct WQuadPolygon
{
	unsigned char* data ;// ËÄÌõ±ßÊ¹ÓÃ1,2,3,4  ÄÚ²¿µã255 Íâ²¿µã0
	int xsize ;
	int ysize ;

	float leftLon , topLat , reso ;
	WQuadPolygon( sllPoint& p0 , sllPoint& p1 , sllPoint& p2 , sllPoint& p3 ,float treso ) ;
	~WQuadPolygon() ;
	void fillLine( sllPoint& p0 , sllPoint& p1 ,unsigned char edgeValue) ;
	
} ;

WQuadPolygon::WQuadPolygon( sllPoint& p0 , sllPoint& p1 , sllPoint& p2 , sllPoint& p3 , float treso )
{
	float lonmin = wmin( wmin( p0.lon , p1.lon ) , wmin(p2.lon,p3.lon ) ) ;
	float latmin = wmin( wmin( p0.lat , p1.lat ) , wmin(p2.lat,p3.lat ) ) ;
	float lonmax = wmax( wmax( p0.lon , p1.lon ) , wmax(p2.lon,p3.lon ) ) ;
	float latmax = wmax( wmax( p0.lat , p1.lat ) , wmax(p2.lat,p3.lat ) ) ;
	this->reso = treso ;

	int tempx0 = lon2grid( lonmin , reso , lonmin - reso/2 ) ;
	int tempx1 = lon2grid( lonmax , reso , lonmin - reso/2 ) ;
	int tempy0 = lat2grid( latmin , reso , latmax + reso/2 ) ;
	int tempy1 = lat2grid( latmax , reso , latmax + reso/2 ) ;

	xsize = tempx1 - tempx0 + 1 ;
	ysize = tempy0 - tempy1 + 1 ;

	leftLon = lonmin ;
	topLat = latmax ;
	
	int asize = xsize * ysize ;
	data = new unsigned char[asize] ;
	memset( data , 0 , asize ) ;

	this->fillLine( p0 , p1,1) ;//edge 1
	this->fillLine( p1 , p2,2) ;//edge 2
	this->fillLine( p2 , p3,3) ;//edge 3
	this->fillLine( p3 , p0,4) ;//edge 4
	
	//Ìî³äÆäËûÎ»ÖÃ
	for(int iy = 0 ; iy < ysize ; ++ iy )
	{
		int vx0 = -1 ;
		int vx1 = -1 ;
		for(int ix = 0 ; ix < xsize ; ++ ix )
		{
			if( data[ix + iy * xsize] )
			{
				if( vx0 < 0 ){
					vx0 = ix ;
				}
				else
				{
					vx1 = ix ;
				}
			}
		}
		for(int ix = vx0+1 ; ix < vx1 ; ++ ix )
		{
			if( data[ix+iy*xsize]==0 )
			{
				data[ix + iy * xsize] = 255 ;
			}
		}
	}

}

WQuadPolygon::~WQuadPolygon()
{
	if( data != 0 )
	{
		delete[] data ;
		data = 0 ;
	}
}
//Ïß¶ÎÀëÉ¢»¯
void WQuadPolygon::fillLine( sllPoint& p0 , sllPoint& p1 ,unsigned char edgeValue ) 
{
	float dlon = p1.lon-p0.lon  ;
	float dlat = p1.lat-p0.lat  ;

	if( dlon == 0 && dlat == 0 )
	{
		int ix = lon2grid( p0.lon ,reso , leftLon - reso/2 );
		int iy = lat2grid( p0.lat , reso , topLat + reso/2 ) ;
		data[ix+iy*xsize] = edgeValue ;
	}
	else if( fabsf(dlon) > fabsf(dlat) )
	{//x·½Ïò
		float k = dlat / dlon  ;
		float b = p0.lat ;
		float clon = p0.lon ;
		float step = reso / 5 ;
		if( dlon < 0 ) step = -step ;
		float tlonmax = p1.lon + step ;
		if( dlon < 0 )
		{
			while( clon > tlonmax )
			{
				float clat = (clon - p0.lon )*k + b ;
				int ix = lon2grid(clon,reso,leftLon - reso/2); //( clon - leftLon ) / reso + 0.5f;
				int iy = lat2grid( clat , reso,topLat + reso/2 );// ( topLat - clat )  / reso + 0.5f;
				data[ix+iy*xsize] = edgeValue ;
				clon += step ;
			}
		}else
		{
			while( clon < tlonmax )
			{
				float clat = (clon - p0.lon )*k + b ;
				int ix =  lon2grid(clon,reso,leftLon - reso/2 );
				int iy = lat2grid( clat , reso,topLat + reso/2 );
				data[ix+iy*xsize] = edgeValue ;
				clon += step ;
			}
		}

	}else
	{
		float k = dlon / dlat ;
		float b = p0.lon ;
		float clat = p0.lat ;
		float step = reso / 5 ;
		if( dlat < 0 ) step = -step ;
		float tlatmax = p1.lat + step ;
		if( dlat < 0 )
		{
			while( clat > tlatmax )
			{
				float clon = (clat - p0.lat )*k + b ;
				int ix =  lon2grid(clon,reso,leftLon -reso/2);
				int iy = lat2grid( clat , reso,topLat +reso/2);
				data[ix+iy*xsize] = edgeValue ;
				clat += step ;
			}
		}else
		{
			while( clat < tlatmax )
			{
				float clon = (clat - p0.lat )*k + b ;
				int ix =  lon2grid(clon,reso,leftLon - reso/2);
				int iy = lat2grid( clat , reso,topLat + reso/2 );
				data[ix+iy*xsize] = edgeValue ;
				clat += step ;
			}
		}
		
	}
}
 

//¼ÆËãÏß¶Î£¨pt0->pt1)ÄÚµÄÏñËØ,½á¹û×·¼Óµ½vecÀïÃæ
//Õâ¸öµØ·½ÐèÒª¿¼ÂÇ¿çÔ½180¶ÈÇé¿ö£¬Á½¸öµã·Ö±æÔÚÍ¼ÏñÁ½²à
void rayline( llPoint& pt0 , llPoint& pt1 , float llreso , 
	float leftlon , float toplat , vector<sllPoint>& vec )
{
	float lon0 = pt0.lon ;
	float lat0 = pt0.lat ;
	float lon1 = pt1.lon ;
	float lat1 = pt1.lat ;

	float dlon = fabsf( lon0-lon1 ) ;
	float dlat = fabsf( lat0-lat1 ) ;

	vec.push_back( sllPoint(lon0 , lat0) ) ;

	int lastix = (lon0 - leftlon) / llreso ;
	int lastiy = (toplat - lat0 ) / llreso ;
	
	if( dlon > dlat )
	{//ÑØ×Åx·½Ïò¼ÆËãÏñËØµã 
		int nstep = fabsf( lon1 - lon0 ) / llreso + 1.f  ;
		float steplon = (lon1-lon0)/nstep ;
		float steplat = (lat1-lat0)/nstep ;
		float clon = lon0 ;
		float clat = lat0 ;
		for(int i = 0 ; i<nstep ; ++ i )
		{
			clon += steplon ;
			clat += steplat ;
			int cx = (clon - leftlon) / llreso  ;
			int cy = (toplat - clat)  / llreso  ;
			if( cx == lastix && cy == lastiy )
			{
				//ÖØ¸´µãÌø¹ý
			}else
			{
				lastix = cx ;
				lastiy = cy ;
				vec.push_back( sllPoint(clon , clat) ) ;
			}
		}

	}else
	{//ÑØ×Åy·½Ïò¼ÆËãÏñËØµã
		int nstep = fabsf( lat1 - lat0 ) / llreso + 1.f ;
		float steplon = (lon1-lon0)/nstep ;
		float steplat = (lat1-lat0)/nstep ;
		float clon = lon0 ;
		float clat = lat0 ;
		for(int i = 0 ; i<nstep ; ++ i )
		{
			clon += steplon ;
			clat += steplat ;
			int cx = (clon - leftlon) / llreso  ;
			int cy = (toplat - clat) / llreso  ;
			if( cx == lastix && cy == lastiy )
			{
				//ÖØ¸´µãÌø¹ý
			}else
			{
				lastix = cx ;
				lastiy = cy ;
				vec.push_back( sllPoint(clon , clat) ) ;
			}
		}
	}

	vec.push_back( sllPoint(lon1 , lat1) ) ;
}


//gridÏñËØµãÅÅÐò
bool wPointCompAsc( const llPoint& left , const llPoint& right )
{
	return  (left.x + left.y * llPoint::gridxsize ) < (right.x + right.y * llPoint::gridxsize ) ;
}


//ÏßÐÔÄÚ²å·µ»Ø xÔÚ (x0,y0) (x1,y1)È·¶¨µÄÖ±ÏßµÄyÖµ
float linear( float x0 , float x1 , float y0 , float y1 , float x )
{
	float sign = (x-x0) * (x-x1) ;
	float y = y0 ;
	if( fabsf(x0-x1) < 0.00001 )
	{
		y = (y0+y1)/2;
	}else
	{
		if( sign > 0 )
		{
			if( x0 >= x1 )
			{
				y = (x-x1)/(x0-x1) * (y0-y1) + y1 ;
			}else
			{
				y = (x-x0)/(x1-x0) * (y1-y0) + y0 ;
			}
		}else
		{
			y = (x-x0)/(x1-x0) * (y1-y0) + y0 ;
		}
	}
	return y ;
}

//Ë«ÏßÐÔ²îÖµ
//pt0 pt1
//pt3 pt2
void bilinear2( float lon , float lat , 
	llPoint& pt0 , 
	llPoint& pt1 , 
	llPoint& pt2 , 
	llPoint& pt3 ,
	vector<float> & res )
{
	static vector<float> inter0 = pt0.values ;
	static vector<float> inter1 = pt0.values ;
	//static vector<float> inter3 = pt0.values ;
	static int nbBil = inter0.size() ;
	float newy0 = 0 ;
	float newy1 = 0 ;

	// interpol lon0 -- lonx -- lon1
	if( fabsf(pt1.lon-pt0.lon) < 0.000001f )
	{//almost zero.
		for(int i = 0 ; i < nbBil ; ++ i )
		{
			inter0[i] = ( pt1.values[i] + pt0.values[i] ) / 2;
		}
		newy0 = (pt1.lat+pt0.lat)/2 ;
	}else
	{
		float dx = ( lon - pt0.lon ) / ( pt1.lon - pt0.lon ) ;
		for(int i = 0 ; i < nbBil ; ++ i )
		{
			inter0[i] = ( pt1.values[i] - pt0.values[i] ) * dx + pt0.values[i] ;
		}
		newy0 = ( pt1.lat - pt0.lat ) * dx + pt0.lat ; 
	}

	if( fabsf(pt3.lon-pt2.lon) < 0.000001f )
	{
		for(int i = 0 ; i < nbBil ; ++ i )
		{
			inter1[i] = ( pt2.values[i] + pt3.values[i] ) / 2;
		}
		newy1 = (pt2.lat+pt3.lat)/2 ;
	}else
	{
		float dx = ( lon - pt3.lon ) / ( pt2.lon - pt3.lon ) ;
		for(int i = 0 ; i < nbBil ; ++ i )
		{
			inter1[i] = ( pt2.values[i] - pt3.values[i] ) * dx + pt3.values[i] ;
		}
		newy1 = ( pt2.lat - pt3.lat ) * dx + pt3.lat ; 
	}

	if( fabsf( newy0 - newy1 )  < 0.000001f )
	{
		for(int i = 0 ; i<nbBil ; ++ i )
		{
			res[i] = ( inter1[i] + inter0[i] ) / 2.f ;
		}
	}else
	{
		float dx = ( lat - newy0 ) / ( newy1 - newy0 ) ;
		for(int i = 0 ; i<nbBil ; ++ i )
		{
			res[i] = ( inter1[i] - inter0[i] ) * dx + inter0[i] ;
		}
	}
}



//µ¥ÏßÐÔ²îÖµ
//pt0 pt1
void singlelinear( float lon , float lat , 
	llPoint& pt0 , 
	llPoint& pt1 ,
	vector<float> & res )
{
	static int nbBil = pt0.values.size() ;

	float dist0 = sqrtf( (lon - pt0.lon )*( lon - pt0.lon ) + (lat - pt0.lat )*(lat - pt0.lat ) ) ;
	float dist1 = sqrtf( (lon - pt1.lon )*( lon - pt1.lon ) + (lat - pt1.lat )*(lat - pt1.lat ) );
	float dista = dist0 + dist1 ;
	if( dista < 0.000001f )
	{//almost zero.
		for(int i = 0 ; i < nbBil ; ++ i )
		{
			res[i] = ( pt1.values[i] + pt0.values[i] ) / 2;
		}
	}else
	{
		float weight0 =  dist0 / dista ;
		float weight1 = 1.f - weight0 ;
		for(int i = 0 ; i < nbBil ; ++ i )
		{
			res[i] = pt1.values[i]*weight1 + pt0.values[i] * weight0;
		}
	}
}


//»ñÈ¡SwathËÄ½Ç¾­Î³¶È·¶Î§
void getSwathLonLatExtent( GDALDataset* dslon , GDALDataset* dslat ,
	double* leftLon , double* rightLon , double* topLat , double* bottomLat )
{
	float lonmin(500) , lonmax(-500) ;
	float latmin(500) , latmax(-500) ;

	int xsize = dslon->GetRasterXSize() ;
	int ysize = dslon->GetRasterYSize() ;

	size_t asize = xsize * ysize ;
	float* lonbuff = new float[asize] ;
	float* latbuff = new float[asize] ;

	dslon->GetRasterBand(1)->RasterIO( GF_Read , 0 , 0 , xsize , ysize , lonbuff , xsize , ysize , GDT_Float32 , 0,0,0) ;
	dslat->GetRasterBand(1)->RasterIO( GF_Read , 0 , 0 , xsize , ysize , latbuff , xsize , ysize , GDT_Float32 , 0,0,0) ;

	for(size_t it = 0 ; it < asize ; ++ it )
	{
		if( lonbuff[it] > -181 && lonbuff[it] <181 
			&& latbuff[it] > -91 && latbuff[it] < 91 )
		{
			if( lonbuff[it] > lonmax ) lonmax = lonbuff[it] ;
			if( lonbuff[it] < lonmin ) lonmin = lonbuff[it] ;
			if( latbuff[it] > latmax ) latmax = latbuff[it] ;
			if( latbuff[it] < latmin ) latmin = latbuff[it] ;
		}

	}

	delete[] lonbuff ;
	delete[] latbuff ;
 
	*leftLon = lonmin  ;
	*rightLon = lonmax ;

	*topLat = latmax ;
	*bottomLat = latmin ;
}
//ÅÐ¶ÏÊÇ·ñ¿ç¹ý180¶È¾­Ïß
bool isCross180Long( llPoint& pt0 , llPoint& pt1 , llPoint& pt2 , llPoint& pt3 )
{
	bool has179 = false ;
	bool has181 = false ;

	if( pt0.lon > 170 && pt0.lon < 181.0 ) has179 = true ;
	else if( pt1.lon > 170 && pt1.lon < 181.0 ) has179 = true ;
	else if( pt2.lon > 170 && pt2.lon < 181.0 ) has179 = true ;
	else if( pt3.lon > 170 && pt3.lon < 181.0 ) has179 = true ;

	if( pt0.lon > -180 && pt0.lon < -170 ) has181 = true ;
	else if( pt1.lon > -181 && pt1.lon < -170 ) has181 = true ;
	else if( pt2.lon > -181 && pt2.lon < -170 ) has181 = true ;
	else if( pt3.lon > -181 && pt3.lon < -170 ) has181 = true ;

	if( has179 && has181 ) return true ;
	else return false ;
}

/// Èç¹û1.4.1²âÊÔÍ¨¹ý£¬ÒÔÏÂ´úÂë¿ÉÒÔÉ¾³ý
//¸ù¾ÝÉý¹ì¡¢½µ¹ì½«¿çÔ½·Ö½çÏßÒÔºóµÄÊý¾ÝÖÃ¿Õ
// ·µ»Ø½á¹û useSide =1Ê¹ÓÃ·Ö½çÏßÉÏÃæ²¿·Ö£¬=-1Ê¹ÓÃ·Ö½çÏßÏÂÃæ²¿·Ö£¬=0 Ã»ÓÐ·Ö½çÏß
// wantAscending = true Ê¹ÓÃÉý¹ì²¿·Ö   =false Ê¹ÓÃ½µ¹ì²¿·Ö
// observeFromUpDown =true ¹ìµÀÊý¾Ý¹Û²âË³Ðò´ÓÉÏ£¨0£©µ½ÏÂ£¨ysize-1£©£¬=false ¹ìµÀÊý¾Ý¹Û²âË³Ðò´ÓÏÂ£¨ysize-1£©µ½ÉÏ£¨0£©
void chooseUseSideForAscendingOrDesending(
	float* latData,
	int* boundaryLineY,
	int xsize,
	int ysize , 
	USETRACKTYPE wantTrack , 
	bool observeFromUpDown , 
	USESIDETYPE* useSide ) 
{
	if( wantTrack == TRACK_BOTH ) return ;
	*useSide= USE_SIDE_FROM_UNKNOWN;
	for(int i = 0 ; i<xsize ; ++ i ) boundaryLineY[i] = -1 ;
	int it = 0;
	if( observeFromUpDown )
	{//¹Û²â·½Ïò´Ó0µ½ysize
		for(int ix = 0 ; ix < xsize ; ++ ix )
		{
			for(int iy = 1 ; iy < ysize -1 ; ++ iy )
			{
				it = ix+iy*xsize ;
				
				if( latData[it] > latData[it-xsize] && latData[it] >= latData[it+xsize] )
				{//×î´óÖµµã
					boundaryLineY[ix] = iy ;
					if( *useSide == 0 )
					{
						if( wantTrack == TRACK_ASC )
						{
							*useSide =  USESIDETYPE::USE_SIDE_FROM_0_TO_ROWX;//1;
						}else
						{
							*useSide = USESIDETYPE::USE_SIDE_FROM_ROWX_TO_YSIZE ;//-1;
						}
					}
					break ;
				}else if( latData[it] < latData[it-xsize] && latData[it] <= latData[it+xsize] ) 
				{//×îÐ¡Öµµã
					boundaryLineY[ix] = iy ;
					if( *useSide == 0 )
					{
						if( wantTrack == TRACK_ASC  )
						{
							*useSide = USESIDETYPE::USE_SIDE_FROM_ROWX_TO_YSIZE ;// -1;
						}else
						{
							*useSide = USESIDETYPE::USE_SIDE_FROM_0_TO_ROWX ;// 1;
						}
					}
					break ;
				}
			}
		}
	}else
	{//¹Û²â·½Ïò´Óysizeµ½0
		for(int ix = 0 ; ix < xsize ; ++ ix )
		{
			for(int iy = 1 ; iy < ysize -1 ; ++ iy )
			{
				it = ix+iy*xsize ;
				if( latData[it] > latData[it-xsize] && latData[it] >= latData[it+xsize] )
				{//×î´óÖµµã
					boundaryLineY[ix] = iy ;
					if( *useSide == 0 )
					{
						if( wantTrack == TRACK_ASC  )
						{
							*useSide = USESIDETYPE::USE_SIDE_FROM_ROWX_TO_YSIZE ;//-1;
						}else
						{
							*useSide = USESIDETYPE::USE_SIDE_FROM_0_TO_ROWX ;// 1;
						}
					}
					break ;
				}else if( latData[it] < latData[it-xsize] && latData[it] <= latData[it+xsize] ) 
				{//×îÐ¡Öµµã
					boundaryLineY[ix] = iy ;
					if( *useSide == 0 )
					{
						if( wantTrack == TRACK_ASC  )
						{
							*useSide = USESIDETYPE::USE_SIDE_FROM_0_TO_ROWX ;//1;
						}else
						{
							*useSide = USESIDETYPE::USE_SIDE_FROM_ROWX_TO_YSIZE ;//-1;
						}
					}
					break ;
				}
			}
		}
	}
}

//ÅÐ¶ÏµãÊÇ·ñÔÚËÄ±ßÐÎÄÚ²¿

sllPoint vect2d(sllPoint& p1, sllPoint& p2) {
    sllPoint temp;
    temp.lon = (p2.lon - p1.lon);
    temp.lat = -1 * (p2.lat - p1.lat);
    return temp;
}
bool pointInRectangle(sllPoint& A, sllPoint& B, sllPoint& C, sllPoint& D, sllPoint& m ) {
    sllPoint AB = vect2d(A, B);  float C1 = -1 * (AB.lat*A.lon + AB.lon*A.lat); float D1 = (AB.lat*m.lon + AB.lon*m.lat) + C1;
    sllPoint AD = vect2d(A, D);  float C2 = -1 * (AD.lat*A.lon + AD.lon*A.lat); float D2 = (AD.lat*m.lon + AD.lon*m.lat) + C2;
    sllPoint BC = vect2d(B, C);  float C3 = -1 * (BC.lat*B.lon + BC.lon*B.lat); float D3 = (BC.lat*m.lon + BC.lon*m.lat) + C3;
    sllPoint CD = vect2d(C, D);  float C4 = -1 * (CD.lat*C.lon + CD.lon*C.lat); float D4 = (CD.lat*m.lon + CD.lon*m.lat) + C4;
    return 0 >= D1 && 0 >= D4 && 0 <= D2 && 0 >= D3;
}

struct IntArray
{
	int* dataPtr ;
	int  arrSize ;
	IntArray( int size ) ;
	inline ~IntArray(){ if( dataPtr ) delete[] dataPtr ; dataPtr = 0 ; } ;
	void fill(int val) ;
private:
	inline IntArray( IntArray& arr ) {};
};
IntArray::IntArray( int size ) 
{
	arrSize = size ;
	dataPtr = new int[arrSize] ;
	if( dataPtr==0 )
	{
		cout<<"IntArray alloc memory failed with size of " << size <<endl ;
		exit(1) ;
	}
}

void IntArray::fill(int val) 
{
	for(int it = 0 ; it < this->arrSize ; ++ it )
	{
		dataPtr[it] = val ;
	}
}


 

string getDirByFilepath (string& str)
{
  size_t found =str.find_last_of("/\\");
  if( found== string::npos) return "" ;
  else
  {
	  return str.substr(0,found+1) ;
  }
}



 
///
int main(int argc, char** argv)
{
	cout<<"A program to convert giirs image into longlat grid image."<<endl ;
	cout<<"wangfengdev@163.com 2020-4-9"<<endl ;
	cout<<"v1.0"<<endl ;//
	cout<<"v1.1 extract scan segments from longitude dataset."<<endl ;//
	cout<<"v1.1.2  use MW long/lat replace LW long/lat 2020-7-9"<<endl ;
	cout<<"v1.2.0 skip 1.1.3,1.1.4 new version support Doc. Lv's Level-3 TSA data.2021-1-13"<<endl ;
	cout<<"v2.0.0 support FY4B GIIRS AVP L2 12km. 2023-4-9"<<endl;
	cout<<"v3.0.1 fea: support FY4B AVP humidity(AQ_Prof) 2024/7/3."<<endl;


	cout<<"syntax:whgiirs2grid input output [dsdata dslon dslat validmin validmax outoffset outscale outfill]"<<endl ;	
	cout<<"syntax:whgiirs2grid fy4b input output"<<endl ;
	cout<<"input 输入文件要求是hdf或者nc格式"<<endl ;
	cout<<"output 输出结果文件（GeoTiff），输出结果为WGS84坐标"<<endl ;
	cout<<"dsdata 数据集在文件中的名称"<<endl ;
	cout<<"dslon  经度数据集在文件中的名称"<<endl ;
	cout<<"dslat  纬度数据集在文件中的名称"<<endl ;
	cout<<"validmin 数据有效区间最小值"<<endl ;
	cout<<"validmax 数据有效区间最大值"<<endl ;
	cout<<"outoffset 输出数据偏移，输出数据公式 outdn=(indata+outoffset)*outscale"<<endl;
	cout<<"outscale 输出数据放大系数"<<endl ;
	cout<<"outfill 输出数据填充值，当输入数据超出有效区间的时候直接使用填充值"<<endl ;
	cout<<"mod1 sample call:whgiirs2grid giirs.hdf output.tif"<<endl ;
	cout<<"mod2 sample call:whgiirs2grid giirsl3.hdf output.tif TSA lon lat -40 40 0 100 -9999"<<endl ;
	cout<<"mod_4b sample call:whgiirs2grid fy4b FY4B...HDF output.tif"<<endl ;
	cout<<"mod_4b_AQ sample call:whgiirs2grid fy4baq FY4B...HDF output.tif"<<endl ;
	cout<<"mod_4b and mod_4b_AQ use 12km as resolution, for FY4A use 16km as resolution."<<endl ;
	cout<<"mod_4b    realvalue=pixelvalue/100.0+100.0 "<<endl;
	cout<<"mod_4b_AQ realvalue=pixelvalue             "<<endl;
	int runMode = MOD1 ;//

	if( argc==3 )
	{
		runMode = MOD1 ;

	}else if( argc == 11 ){
		//good
		runMode = MOD2 ;
	}
	else if( argc==4 )
	{
		if( strcmp(argv[1],"fy4b")==0 ) {
			runMode = MOD_4B;
		}else if(strcmp(argv[1],"fy4baq")==0) {
			runMode = MOD_4B_AQ;
		}	
		else{
			cout<<"unsupported mode: "<< argv[1] <<endl ;
			return 11 ;
		}
		
	}
	else
	{
		cout<<"params should be 3, 4 or 11."<<endl ;
		return 11 ;
	}

	cout<<"run mode:"<<runMode<<endl ;

	string dsName = "AT_Prof" ;
	string lonName = "MW_Longitude" ;
	string latName = "MW_Latitude" ;

	double inValidMin = 0 ;
	double inValidMax = 500 ;
	double outOffset = -100 ;
	double outScale = 100 ;
	double outFill = -9999 ;

	string inHdfFile = "";
	string outfile = "" ;

	if( runMode==MOD1 )
	{
		inHdfFile = argv[1] ;
		outfile = argv[2] ;
	}
	if( runMode==MOD2 ){

		inHdfFile=argv[1] ;
		outfile = argv[2] ;

		dsName = argv[3] ;
		lonName = argv[4] ;
		latName = argv[5] ;

		inValidMin = atof( argv[6] ) ;
		inValidMax = atof( argv[7] ) ;
		outOffset = atof( argv[8] ) ;
		outScale = atof( argv[9] ) ;
		outFill = atof( argv[10] ) ;
	}else if( runMode==MOD_4B ){
		dsName = "AT_Prof" ;
		inHdfFile = argv[2] ;
		outfile = argv[3] ;
		outOffset = -100 ;
		outScale = 100 ;
	}else if( runMode==MOD_4B_AQ ) {
		dsName = "AQ_Prof" ;
		inHdfFile = argv[2] ;
		outfile = argv[3] ;
		outOffset = 0 ;
		outScale =  1.0 ;
	}

	cout<<"inHdfFile:"<<inHdfFile<<endl ;
	cout<<"outfile:"<<outfile<<endl;

	cout<<"info dsName:"<<dsName<<endl ;
	cout<<"info lonName:"<<lonName<<endl ;
	cout<<"info latName:"<<latName<<endl ;
	cout<<"info inValidMin:"<<inValidMin<<endl ;
	cout<<"info inValidMax:"<<inValidMax<<endl ;
	cout<<"info outOffset:"<<outOffset<<endl ;
	cout<<"info outScale:"<<outScale<<endl ;
	cout<<"info outFill:"<<outFill<<endl ;



	string datafile = string("HDF5:") + inHdfFile + "://" + dsName ;
	string lonfile = string("HDF5:") +  inHdfFile + "://" + lonName ;//bugfixed use MW replace LW 2020-7-9
	string latfile = string("HDF5:") +  inHdfFile + "://" + latName ;//bugfixed use MW replace LW 2020-7-9
 

	double reso = 0.16 ;//FY4A
	double leftLon = 66.0;
	double rightLon = 144.0;
	double topLat = 55.0;
	double bottomLat = 3.0;

	if( runMode==MOD_4B || runMode==MOD_4B_AQ){
		reso = 0.12 ;//FY4B
		leftLon = 66.0;
		rightLon = 144.0;
		topLat = 55.0;
		bottomLat = 3.0;
	}
	cout<<"output resolution:"<<reso<<endl ;
	cout<<"output left,right top,bottom:"<<leftLon<<","<<rightLon<<" "<<topLat<<","<<bottomLat<<endl;


	// double inFillValue = -999999 ;//ÓÐÐ§Öµ·¶Î§0-500 float32
	// double outFillValue = -9999 ;//Êä³öÖµ 273.51 -> 273.51*100 -> 27351 -10000 -> 17351 short

	GDALAllRegister() ;

	wGdalRasterFloat rasterlon , rasterlat , rasterdata ;
	bool oklon = rasterlon.open( lonfile ) ;
	bool oklat = rasterlat.open( latfile ) ;
	bool okdata = rasterdata.open( datafile ) ;

	if( oklon == false )
	{
		cout<<"failed to open "<<lonfile<<endl ;
		return 12 ;
	}
	if( oklat == false )
	{
		cout<<"failed to open "<<latfile<<endl ;
		return 13 ;
	}
	if( okdata == false )
	{
		cout<<"failed to open "<< datafile<<endl ;
		return 14 ;
	}
  

	int nband = rasterdata.getNBand() ;
	size_t datacol = rasterdata.getXSize() ;
	size_t datarow = rasterdata.getYSize() ;
	const int ncxsize = datacol ;
	const int ncysize = datarow ; 

	cout<<"data-size ncol nrow nband:"<<datacol<<" "<<datarow<<" "<<nband<<endl ;


	int leftLonix =lon2grid(leftLon,reso , -180 ) ;
	int rightLonix =lon2grid(rightLon , reso , -180 );

	leftLon = grid2lon(leftLonix,reso , -180) ;
	rightLon = grid2lon(rightLonix,reso , -180) ;

	int topLatiy = lat2grid( topLat , reso , 90 ) ;
	int bottomLatiy = lat2grid( bottomLat , reso , 90) ;

	topLat = grid2lat( topLatiy , reso , 90) ;
	bottomLat = grid2lat( bottomLatiy , reso , 90) ;

	size_t gridXSize = rightLonix - leftLonix + 1 ;
	size_t gridYSize = bottomLatiy - topLatiy + 1 ;
	llPoint::gridxsize = gridXSize ;

	wGdalRasterI16 rasterOut ;
	wGdalRasterFloat rasterOut3 ;//for MOD_4B_AQ
	bool outok = rasterOut.create( gridXSize , gridYSize , nband ) ;
	if( outok==false )
	{
		cout<<"failed to create rasterOut."<<endl ;
		return 14 ;
	}
	outok = rasterOut3.create( gridXSize , gridYSize , nband ) ;


	for(int i = 0 ; i<nband ; ++ i )
	{
		rasterOut.fill(i,outFill) ;
		rasterOut3.fill( i, outFill) ;
	}

	double trans[6] ;
	trans[0] = leftLon ;
	trans[1] = reso ;
	trans[2] = 0 ;
	trans[3] = topLat ;
	trans[4] = 0 ;
	trans[5] = -reso ;
	rasterOut.copyTrans( trans ) ;
	rasterOut.copyProj("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,"
		"AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,"
		"AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //wgs84

	rasterOut3.copyTrans( trans ) ;
	rasterOut3.copyProj("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,"
		"AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,"
		"AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //wgs84



	int per0 = -1 ;
	cout<<"processing..."<<endl ;
	
	vector<int> scanStartColVector ;
	scanStartColVector.push_back(0) ;

	if( MOD1==runMode || MOD2==runMode )
	{//For FY4A
		int lastrowx0 = 31*ncxsize ;//FY4A
		for(int ix = 0 ; ix < ncxsize-1 ;++ ix )
		{
			float lon1 = rasterlon.getValuef( ix + lastrowx0,0) ;
			float lon2 = rasterlon.getValuef( ix + lastrowx0 + 1,0) ;
			if( lon1 > -180 && lon1 < 180 )
			{//valid lon
				if( lon2 > -180 && lon2 < 180  )
				{
					if( lon1 > lon2 )
					{
						scanStartColVector.push_back(ix+1) ;
					}
				}else
				{
					scanStartColVector.push_back(ix+1) ;
				}
			}
		}
	}


	cout<<"print segments:"<<endl ;
	for(int iscan = 0 ; iscan <  scanStartColVector.size()  ;++iscan )
	{
		cout<<scanStartColVector[iscan]<<endl ;
	}

	int numScan = scanStartColVector.size()  ;
	int allrows = numScan*32 ;//FY4A
	int rowsPerScan = 32 ;//FY4A

	if( runMode==MOD_4B || runMode==MOD_4B_AQ ){
		allrows = datarow ;//FY4B
		rowsPerScan = datarow ;
	}

	int processRows = 0 ;
	// const int val2outScale = 100 ;
	// const int val2outOffset = -10000 ; inval0*val2outScale + val2outOffset
	for( int iscan = 0 ; iscan < numScan ; ++ iscan )
	{
		int icolstart = scanStartColVector[iscan] ;
		int icolend = 0 ;
		if( iscan < numScan-1 )
		{
			icolend = scanStartColVector[iscan+1] ;
		}else{
			icolend = ncxsize ;
		}
		cout<<"scan:"<<iscan<<" startcol:"<<icolstart<<" endcol:"<<icolend<<endl ;
		for(int irow = 0 ; irow < rowsPerScan-1 ; ++ irow )
		{
			++ processRows ;
			for(int icol = icolstart ; icol < icolend-1 ; ++ icol )
			{
				
				float currLeftLon = leftLon ;
				llPoint pt0 , pt1 , pt2 , pt3 ;
				// pt0 pt1
				// pt3 pt2
				int it = icol + irow * datacol ;
				//ÅÐ¶ÏÊÇ·ñÓÐÎÞÐ§µÄ¾­Î³¶È£¬Èç¹ûËÄ¸öµãÓÐÎÞÐ§µÄ¾­Î³¶ÈÖ±½ÓÌø¹ý
				if( rasterlon.getValuef(it , 0) < leftLon           || rasterlon.getValuef(it,0) > rightLon ) continue ;
				if( rasterlon.getValuef(it+1 , 0) < leftLon         || rasterlon.getValuef(it+1,0) > rightLon) continue ;
				if( rasterlon.getValuef(it+1+datacol , 0) < leftLon || rasterlon.getValuef(it+1+datacol,0) > rightLon) continue ;
				if( rasterlon.getValuef(it+datacol , 0) < leftLon   || rasterlon.getValuef(it+datacol,0) > rightLon) continue ;

				if( rasterlat.getValuef(it , 0) < bottomLat           || rasterlat.getValuef(it,0)              > topLat ) continue ;
				if( rasterlat.getValuef(it+1 , 0) < bottomLat         || rasterlat.getValuef(it+1,0)            > topLat) continue ;
				if( rasterlat.getValuef(it+1+datacol , 0) < bottomLat || rasterlat.getValuef(it+1+datacol,0)    > topLat) continue ;
				if( rasterlat.getValuef(it+datacol , 0) < bottomLat    || rasterlat.getValuef(it+datacol,0)     > topLat) continue ;

				pt0.lon = lon2gridlon( rasterlon.getValuef(it , 0) ,reso) ;
				pt0.lat = lat2gridlat( rasterlat.getValuef(it , 0),reso) ;
				pt0.x = lon2grid(pt0.lon , reso , currLeftLon - reso/2 ) ; //( pt0.lon - leftLon ) / reso  ;
				pt0.y = lat2grid(pt0.lat , reso , topLat + reso/2 ) ;// ( topLat - pt0.lat ) /reso  ;

				pt1.lon = lon2gridlon( rasterlon.getValuef(it+1,0) ,reso );
				pt1.lat = lat2gridlat( rasterlat.getValuef(it+1,0),reso) ;
				pt1.x =  lon2grid(pt1.lon , reso, currLeftLon- reso/2)  ;
				pt1.y =  lat2grid(pt1.lat , reso, topLat+ reso/2)  ;

				pt2.lon =  lon2gridlon( rasterlon.getValuef(it+datacol+1,0),reso) ;
				pt2.lat = lat2gridlat( rasterlat.getValuef(it+datacol+1,0),reso) ;
				pt2.x =  lon2grid(pt2.lon , reso, currLeftLon- reso/2)  ;
				pt2.y =  lat2grid(pt2.lat , reso, topLat+ reso/2)    ;

				pt3.lon =  lon2gridlon( rasterlon.getValuef(it+datacol,0),reso) ;
				pt3.lat =lat2gridlat(   rasterlat.getValuef(it+datacol,0),reso) ;
				pt3.x = lon2grid(pt3.lon , reso, currLeftLon- reso/2)   ;
				pt3.y = lat2grid(pt3.lat , reso, topLat+ reso/2)   ;

				sllPoint sp0(pt0.lon , pt0.lat);
				sllPoint sp1(pt1.lon , pt1.lat);
				sllPoint sp2(pt2.lon , pt2.lat);
				sllPoint sp3(pt3.lon , pt3.lat);
				WQuadPolygon wpoly( sp0 , sp1 , sp2 , sp3 , reso ) ;
				int quadLeft = lon2grid( wpoly.leftLon , reso , currLeftLon - reso/2 ) ;// ( wpoly.leftLon - currLeftLon ) / reso ;
				int quadTop =  lat2grid( wpoly.topLat , reso , topLat + reso/2 );//(topLat - wpoly.topLat ) / reso ;

				//°ÑÊäÈëÖµÐ´Èë½á¹ûºÍ²îÖµµãÖÐ
				for(int ib = 0 ; ib < nband ; ++ ib )
				{
					bool hasFillValue = false ;//Èç¹ûÓÐÌî³äÖµÄÇÃ´²»½øÐÐ²åÖµ¼ÆËã
					
					//µÚÒ»±éÖ±½ÓÐ´ÈëÖµ£¬µÚ¶þ±éÄÚ²å£¬µÚÈý±éÌî²¹ËÄÖÜÏñËØÓÐÖµµÄ¿Õµã
					float inval0 = rasterdata.getValuef(it , ib) ; 
					float inval1 = rasterdata.getValuef(it+1 , ib) ; 
					float inval2 = rasterdata.getValuef(it+datacol+1 , ib) ; 
					float inval3 = rasterdata.getValuef(it+datacol , ib) ; 

					if( inval0 > inValidMin && inval0 < inValidMax )
					{
						rasterOut.setValuei( pt0.x , pt0.y , ib , (inval0+outOffset)*outScale ) ;
						rasterOut3.setValuef( pt0.x , pt0.y , ib , (inval0+outOffset)*outScale ) ;
					}
					else
					{
						hasFillValue = true ;
						rasterOut.setValuei( pt0.x , pt0.y , ib , outFill ) ;
						rasterOut3.setValuef( pt0.x , pt0.y , ib , outFill ) ;
					}

					if( inval1 > inValidMin && inval1 < inValidMax )
					{
						rasterOut.setValuei( pt1.x , pt1.y , ib , (inval0+outOffset)*outScale ) ;
						rasterOut3.setValuef( pt1.x , pt1.y , ib , (inval0+outOffset)*outScale ) ;
					}
					else
					{
						hasFillValue = true ;
						rasterOut.setValuei( pt1.x , pt1.y , ib , outFill ) ;
						rasterOut3.setValuef( pt1.x , pt1.y , ib , outFill ) ;
					}

					if( inval2 > inValidMin && inval2 < inValidMax )
					{
						rasterOut.setValuei( pt2.x , pt2.y , ib , (inval0+outOffset)*outScale ) ;
						rasterOut3.setValuef( pt2.x , pt2.y , ib , (inval0+outOffset)*outScale ) ;
					}
					else
					{
						hasFillValue = true ;
						rasterOut.setValuei( pt2.x , pt2.y , ib , outFill ) ;
						rasterOut3.setValuef( pt2.x , pt2.y , ib , outFill ) ;
					}

					if( inval3 > inValidMin && inval3 < inValidMax )
					{
						rasterOut.setValuei( pt3.x , pt3.y , ib , (inval0+outOffset)*outScale ) ;
						rasterOut3.setValuef( pt3.x , pt3.y , ib , (inval0+outOffset)*outScale ) ;
					}
					else
					{
						hasFillValue = true ;
						rasterOut.setValuei( pt3.x , pt3.y , ib , outFill ) ;
						rasterOut3.setValuef( pt3.x , pt3.y , ib , outFill ) ;
					}
					//pt0.values.push_back( theDataBuffer.getScanValue(ib,it) );
					if( hasFillValue ) 
					{
						continue ;//ÓÐÎÞÐ§Öµ²»×öÄÚ²å£¬½øÐÐÏÂÒ»²¨¶Î¼ÆËã
					}
					pt0.values.clear() ;
					pt1.values.clear() ;
					pt2.values.clear() ;
					pt3.values.clear() ;
					pt0.values.push_back( inval0 );
					pt1.values.push_back( inval1 );
					pt2.values.push_back( inval2 );
					pt3.values.push_back( inval3 );

					//¿ªÊ¼ÄÚ²å
					for(int iby = 0 ; iby < wpoly.ysize ; ++ iby )
					{
						for(int ibx = 0 ; ibx < wpoly.xsize ; ++ ibx )
						{
							int quadPixelType = wpoly.data[ibx + iby * wpoly.xsize ] ;
							if( quadPixelType >0 )
							{//ÔÚËÄ±ßÐÎÄÚ²¿
								vector<float> bilresult = pt0.values ;
								double lon1 = grid2lon( ibx , reso , wpoly.leftLon- reso/2) ;
								double lat1 = grid2lat( iby , reso , wpoly.topLat+ reso/2 ) ;
								int gridPtx = quadLeft + ibx  ;//( lon1 - currLeftLon ) / reso  ;
								int gridPty = quadTop + iby ;//( topLat - lat1 )  / reso  ;

								int tempgridit0 = gridPty * gridXSize + gridPtx ;
								if( rasterOut.getValuei(tempgridit0 , ib) != outFill ) continue ; //¸ÃÏñËØÒÑ¾­ÓÐÖµÁË£¬Ìø¹ý£¬¼ÆËãÏÂÒ»¸ö
								if( rasterOut3.getValuei(tempgridit0 , ib) != outFill ) continue ;
							
								if( quadPixelType == 255 )
								{
									bilinear2( lon1  , lat1 ,
										pt0,pt1,pt2,pt3 , bilresult
										) ;
								
								}else if( quadPixelType == 1 )
								{//±ßÔµÏßÉÏÒ»Î¬¾àÀë²åÖµ
									singlelinear( lon1 , lat1 , pt0 , pt1 , bilresult ) ;
								}else if( quadPixelType == 2 )
								{
									singlelinear( lon1 , lat1 , pt1 , pt2 , bilresult ) ;
								}else if( quadPixelType == 3 )
								{
									singlelinear( lon1 , lat1 , pt2 , pt3 , bilresult ) ;
								}else// if( wpoly.data[ibx + iby * wpoly.xsize ] == 4 )
								{
									singlelinear( lon1 , lat1 , pt3 , pt0 , bilresult ) ;
								} 
								rasterOut.setValuei( tempgridit0 , ib , (bilresult[0]+outOffset)*outScale ) ;
								rasterOut3.setValuef( tempgridit0 , ib , (bilresult[0]+outOffset)*outScale ) ;
							}
						}
					}//end for iby
				}//end for ib
			}//end icol
			int per1 = processRows * 100 / allrows ;
			if( per0 != per1 )
			{
				per0 = per1 ;
				cout<<per1<<" " ;
			}
		}//end irow
	}
	cout<<"100."<<endl ;
	per0 = -1 ;
	cout<<"fill gaps."<<endl ;
	for(size_t irow = 1 ; irow <  gridYSize-1 ; ++ irow )
	{
		for( size_t icol = 1 ; icol < gridXSize-1 ; ++ icol )
		{
			size_t it = irow * gridXSize + icol ;
			for(int ib = 0 ; ib < nband ; ++ ib )
			{
				if( rasterOut.getValuei(it , ib ) == outFill )
				{
					int upval = rasterOut.getValuei(it-gridXSize , ib )  ;
					int downval = rasterOut.getValuei(it+gridXSize , ib ) ;
					int leftval = rasterOut.getValuei(it-1 , ib ) ;
					int rightval = rasterOut.getValuei(it+1 , ib ) ;
					if( upval != outFill && downval != outFill && leftval != outFill && rightval != outFill )
					{
						int averval = (upval + downval + leftval + rightval)/4 ;
						rasterOut.setValuei(it , ib , averval ) ;
					}
				}
				if( rasterOut3.getValuei(it , ib ) == outFill )
				{
					float upval = rasterOut3.getValuef(it-gridXSize , ib )  ;
					float downval = rasterOut3.getValuef(it+gridXSize , ib ) ;
					float leftval = rasterOut3.getValuef(it-1 , ib ) ;
					float rightval = rasterOut3.getValuef(it+1 , ib ) ;
					if( (int)upval != outFill && (int)downval != outFill && (int)leftval != outFill && (int)rightval != outFill )
					{
						float averval = (upval + downval + leftval + rightval)/4 ;
						rasterOut3.setValuef(it , ib , averval ) ;
					}
				}
			}
			
		}
		int per1 = irow * 100.f / gridYSize ;
		if( per1!=per0 )
		{
			per0 = per1 ;
			cout<<per0<<"." ;
		}
	}
	cout<<"100."<<endl ;

	cout<<"fill gaps end."<<endl ;
	if( runMode==MOD_4B_AQ ) {
		rasterOut3.save( outfile ) ;
	}else{
		rasterOut.save( outfile ) ;
	}

	return 0;
}

